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ABSTRACT 

The presence of a second planet in a known, transiting-planet system will cause the 
time between transits to vary. These variations can be used to constrain the orbital 
elements and mass of th e perturbing planet. We an alyse the set of transit times of the 
TrES-1 system given in ICharbonneau et alJ (|2005f) . We find no convincing evidence 
for a second planet in the TrES-1 system from that data. By further analysis, we 
constrain the mass that a perturbing planet could have as a function of the semi- 
major axis ratio of the two planets and the eccentricity of the perturbing planet. Near 
low-order, mean-motion resonances (within ~ 1% fractional deviation), we find that 
a secondary planet must generally have a mass comparable to or less than the mass 
of the Earth-showing that this data is the first to have sensitivity to sub Earth-mass 
planets. We compare the sensitivity of this technique to the mass of the perturbing 
planet with future, high-precision radial velocity measurements. 
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1 INTRODUCTION 

Since the first discovery, a decade ago, of a planet orbitin g 
a distant main- sequence star (iMavor fc Quelozl Il995l) . 
several search techniques have been employed to identify 
additional extrasolar planets. These techniqu es include 



radial vel ocity (iMarcv et all 200? ; Moutou et"al"ll200 



trometry (|FoiTLfc_Jkernato 



2003 



Glindemann et al 



XL 
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2003; 



iPravdo et al jl2004l : iGouda et all2005t|Lattanzi et al.l2005l) 
plane tary microlensing - ^U^r^^^^l^^^^]^^^^^^^L| 

20031 ). and p lane t ary transits ([K pnacki et alj .2003: 



Bou chy et 
2005j; |j 



id plane 
all 12004 



Konacki et al 



Pont et afl 12004 ICharbonneau et all 
A recent addition to the 
repertoire of planet search techniques consists of looking 
for additional planets in a known, transiting system by 
analyzing the variation in the time between planetary 
transits. These transit timing variations (TTV) can be 
used to constrain the orbital elements of an unseen, 
perturbing planet, even if its mass is comparable to the 
mass of the Earth jMiralda-Escudel 120021: |Agol et alJl2005t 
iHolman fc Murravl2005h . Thus, for the near term, TTV can 
detect planets around main-sequence stars that are too small 
to detect by any other means. This sensitivity is particularly 
manifest near mean-motion resonances, which recent works 



by iThommesl i2005l) . iPapaloizou fc Szuszkiewiczl (2005 
and lZhou et alJ 720051) suggest might be very common. 

The study of extrasolar planetary systems allows us to 
address fundamental questions about mechanisms of planet 
formation, the prevalence of small, rocky planets like our 
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Earth, and the evolution of multiple planetary systems. Mul- 
tiple planet systems are an important subset of extrasolar 
planets because they give insight into our own solar system. 
TTV should prove to be an important tool in sampling the 
population of such planets since it has the ability to discern 
the presence of bodies with very little mass (< M©). 



The application of TTV depends upon the dis- 
covery and monitoring of transiting planetary sys- 
tems. The fir st successful dete c tion of a plane- 
tary transit iCharbonneau et alJ l200d: iHenrv et alJ 
2000) was for a plane t iden tified spectroscopically, 
HD209458b jMazeh et alJ l2000l) . The first extrasolar 
planet to be discovered from t ransit data was th e 
OGLE-TR-56b system reported bv iKonacki et al] J20031) . 
Existing and future searches for planetary t r ansits , such 
as the COROT faourde. Rouan. fc Leeerl 120031). XO 
jMcCulloueh et al.l 120051) . and Kepler feorucki and et al] 
l20u^ir "missions. are expected to provide many candidate 
systems where the timing of the transits can be analysed. 



The TrES-1 planetary system, reported bv lAlonso et al] 

(2004), was al so discovered via planetary transits. In a re- 
cent paper bv ICharbonneau et al] B2005i) (hereafter C05), 
the detection of thermal emission from the surface of TrES- 
lb was announced. Table 1 of C05 gives the timing of 12 
planetary transits. We analyse the deviations from periodic- 
ity in those data in order to identify and constrain potential 
companion planets in that system. 
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2 THE DATA AND SOFTWARE 



The timing data re ported in C05 were derived from the 11 
transits reported bv lAlonso et al.l (120041) with an additional 
transit that was ob served at the IAC 80cm telescope after 
lAlonso et alJ i2004h went to press. One transit was excluded 
from their analysis because it constituted a 6-a departure 
from a constant period and because of anomalous features 
in the ingress and egress. That point, if it is valid, is the 
most interesting point for our purposes because the TTV 
signal is defined by such deviations. Consequently, we anal- 
yse two different sets of data from C05; the "12-point" set 
which includes this point, and the "11-point" set which does 
not. This s t udy is the f irst analysis of TTV as pre sented by 
lAeol et alJ J2005I) and iHolman fc Murravl J2005I) . And, it 
may be the first search for planets around main-sequence 
stars that can probe masses smaller than the mass of the 
Earth. 

To calculate the time of transit, we integrate the equa- 
tions of motion using the Burlisch-Stoer integration algo- 
rithm. When the transiting planet crosses the line of sight, 
the transit time is identified and tabulated. These simulated 
transit times are compared to the observed transit times and 
the quality of the fit of the orbital elements of the simulated 
system with those of the actual system is characterized by 
the x 2 statistic. For this work we restrict our analysis to 
planets whose orbits are coplanar and edge-on. 

To locate the set of orbital elements that gives the mini- 
mum x 2 , we sequentially employ a random and then a direct 
minimization algorithm. For the random portion we gener- 
ate many sets of orbital elements and record the \ 2 f° r each. 
All parameters are drawn from a uniform distribution, ex- 
cept for the eccentricity and mass which are logarithmically 
uniform; between 10 -3 ' 4 and 1.0 for the eccentricity and be- 
tween 10 -6 and 10 -2 ' 3 for the ratio of the planet to stellar 
mass. We chose the logarithmic distribution so as to not as- 
sume a mass scale a priori-ova limits corresponding approx- 
imately to masses too small to be seen and too large to have 
not been seen. This approach agree s with the empirical dis- 
tribut ion of masses determined by iTabachnik fc Tremaind 
(2002). The logarithmic distribution was chosen for eccen- 
tricity because tidal circularization and orbital stability con- 
straints favor small eccentricities (< 0.1). For completeness, 
we also conducted a study where the eccentricity of the per- 
turbing planet began at zero. 

We tabulate the two systems with the lowest x 2 as start- 
ing points for the direct minimization algorithm. Experi- 
ments with simulated data showed no significant improve- 
ment in the final x 2 when more than two solutions were 
retained. We assume, following the random minimization, 
that at least one set of orbital elements is near the global 
minimum of the x 2 . That is, the nearest local minimum is 
the global minimum and that the local gradient points to- 
ward that minimum. Tests of this approach with simulated 
data confirm this assumption. From these points we use a 
direct minimization algorithm to locate the bottom of the 
global minimum. 



3 RESULTS 

3.1 Search for Secondary Planets 

We conducted a variety of searches for the best-fitting per- 
turbing planet in the TrES-1 system. These searches in- 
cluded different combinations of orbital elements for TrES- 
lb. We found that any reduction in the overall \ 2 obtained 
by including the parameters e and w for TrES-lb was offset 
by the loss of a degree of freedom. Therefore, we report re- 
sults from the search where the eccentricity of TrES-lb was 
fixed at zero. 

We stepped through the semi-major axis ratio of the 
putative secondary planet and TrES-lb. At each point we 
minimized over six parameters: the eccentricity, longitude 
of pericentre, time of pericentre passage, and mass of the 
secondary planet and the period and the initial longitude 
of TrES-lb. The inclination and ascending node of the per- 
turbing planet were identical to the values for TrES-lb. We 
analysed the data for both interior and exterior perturbers. 

Our analysis did not produce any promising solution; 
though we present one interesting case for an interior per- 
turber found from the 12-point analysis. This solution is 
not near a low-order mean-motion resonance. Indeed, the 
X 2 was generally much higher near the j:j + 1 and the j:l 
resonances than in the regions between them. Figure^com- 
pares the timing residuals for this solution with the data. 
The reduced x 2 f° r this system is 2.8 on N — 6 degrees of 
freedom (where N = 12 is the number of data) compared 
with 6.3 for no perturber (N — 2). 

This solution, while it would have been detected from 
RV measurements, is interesting because the average size 
of the timing deviations is larger than that of the data- 
making the variations easier to detect. However, we suspect 
that while this solution is numerically valid, it is merely an 
artifact of the gap between the two primary epochs of ob- 
servation. For this solution, and others like it, the simulated 
timing residuals consist of small, short-term variations su- 
perposed on a large-amplitude, long-period variation with a 
period that is a multiple of the difference between the two 
epochs. Several candiate systems for both the 11-point and 
the 12-point analyses showed similar behaviour. Additional 
data, taken at a time that is not commensurate with this 
period would remove false solutions of this type. 

From the above analysis, we conclude that there is not 
sufficient information in the data to uniquely and satisfac- 
torally determine the characteristics of a secondary planet 
in the TrES-1 system. This is in part because the number 
of model parameters is not much larger than the number 
of data and because the typical timing precision, ~ 100s, is 
not a sufficiently small fraction of the orbital period of the 
transiting planet (about 4 x 10~ 4 ) to distinguish between 
different solutions. In addition, the gap in coverage appears 
to affect the minimization dramatically. We believe that it is 
primarily responsible for the observed fact that two nearby 
sets of orbital elements will have very different values of 
X 2 -a slight change in the long-term variation will cause the 
simulated transit times to miss the second epoch of observa- 
tions. Additional timing data, with precision comparable to 
the most precise of the given data, ~ 30s, and at an epoch 
that is not commensurate with the existing gap in coverage 
will allow for a more complete investigation of the system. 
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Figure 1. A comparison of the TrES-lb transit timing residuals 
to those from the 12-point analysis for an interior perturber. The 
diamonds are the timing residual at each transit, the filled dia- 
monds denote transits where there is data, the squares with error 
bars are the data. The f 2-point solution has a reduced \ 2 of 2.8 on 
6 degrees of freedom with ephemeris parameters m = O.llfMj, 
P = 1.76(2, e = 0.153, w = 154.8°, r = 2452847.4392d. 



3.2 Constraints on Secondary Planets 

That the results of the planet search were inconclusive is not 
particularly surprising since the transit timings show little 
variation from a constant period and the point that deviates 
the most is suspect. If a satisfactory solution cannot be de- 
termined from these data, we can still use them to place con- 
straints on the orbital elements of a two-planet system. Of 
particular interest is a constraint on the mass that the sec- 
ondary planet could have as a function of various orbital ele- 
ments; similar to the constraint that resulted f rom analyzing 
the O GLE-1998-BUL-14 microlensing event (lAlbrow et alJ 
2000). As the mass of the hypothetical perturbing planet 
increases, the x 2 OI t ne model should grow significantly re- 
gardless of the values of the remaining orbital elements. 
Therefore, we systematically studied a grid of values on the 
semi-major axis/eccentricity plane of the perturbing planet- 
assuming that the orbit of TrES-lb is initially circular. At 
each point we identified the mass that a perturber needs in 
order to produce a large deviation from the data. 

At each location in the a/e plane we set the mass of 
the companion to be very small (<C Mg) and, for a random 
value for the longitude of pericenter and the time of pericen- 
ter passage, we calculated the \ 2 of the timing residuals. We 
increased the mass of the secondary planet until the x 2 grew 
by some fiducial amount. At that point we minimized the 
X 2 over the longitude of pericentre and the time of pericen- 
tre passage of the secondary planet and the initial longitude 
of TrES-lb. Following the minimization the x 2 typically fell 
below the fiducial amount and the mass of the perturber 
was again increased. However, if after the minimization the 
X 2 remained above the threshold, the orbital elements of the 
system were recorded, the x 2 threshold increased, and the 
process repeated until either a maximum mass of the per- 
turbing planet or a maximum x 2 was achieved. This pro- 
cedure gives the minimum x 2 as a function of perturbing 
planet mass for each point in the a/e plane. 

The mass that corresponds to a 3-a increase in the 
X 2 constitutes our estimate for the maximum allowed mass 
of a secondary pl anet. We use the approach outlined by 
IPress et alJ dl992t) where we locate the mass that causes 



the difference between the x 2 of the null hypothesis and the 
X 2 obtained with a secondary planet to equal Ax 2 = 9. This 
gives the 3-ct limit because we allow only one parameter to 
vary-the remaining parameters are either fixed or marginal- 
ized at each point. This maximum mass is shown as a func- 
tion of a and e for the 11-point set in Figure for both an 
interior and an exterior perturber. The contours correspond 
to 100, 10, and 1 M®. The dark portion in the upper-left 
corners are where the orbits overlap and we assign a mass of 
1O _8 M0 to those locations. We see from these figures many 
regions where the mass of a companion must be compara- 
ble to or less than the mass of the Earth regardless of its 
orbital eccentricity. The most stringent constraints are near 
the j:j + 1 mean-motion resonances. Of particular interest 
are the very tight constraints for low-eccentricity perturbers 
(which constraints continue to an eccentricity of zero) be- 
cause tidal circularization drives the eccentricity down. The 
12-point analysis gave similar results which are not shown. 



3.3 Comparison with Radial Velocity 

Much of the region where a secondary planet is not tightly 
constrained by our analysis could be limited more efficiently 
by radial velocity measurements. Figure |3] shows the con- 
straint achieved from this TTV analysis and the constraint 
from a hypothetical RV analysis, with the same number of 
data, as a function of the period ratio of the two planets. 
We assume that the RV measurements have a precision of 
both 5 and 1 m s _1 and that the orbit of the perturber can 
be treated as circular. Figures 2] and |S] are the same com- 
parison focussed on the region surrounding the 2:1 and 3:2 
mean-motion resonances respectively. 

We derived the RV limit in these figures by assuming 
that the velocity residuals that remain after removing the ef- 
fects of TrES-lb randomly sample the phase of the putative 
second planet and that the noise is uncorrelated with the 
RV measurements. If no secondary planet exists, then the 
residuals should surround zero with a variance equal to a 2 
where a is the precision of the measurements. The expected 
X 2 is then equal to Xo — N ~ Qo where Qo is the number of 
model parameters. On the other hand, if a secondary planet 
exists, the expected x 2 is 



2 (N-Q) 



2 / 2nGM ^ 2/3 



(1) 



where jj, is the planet to star mass ratio, Mo is the mass of 
the star, G is Newton's constant, and Q is the number of 
model parameters. By subtracting Xo from Xp and solving 
for the mass ratio we obtain the mass of a secondary planet 
that is detectable with a given Ax 2 threshold as a function 
of the period of the planet and the precision of the radial 
velocity measurements 



2(A X 2 + Q- Qo) 



-, 1/2 



N- 



2-kGMo 



1/3 



(2) 



We see from Figure that for low-order, mean-motion 
resonances the TTV analysis is able to place constraints on 
the mass that are nearly an order-of-magnitude smaller than 
the RV technique with 5 m s" 1 precision and that there are 
regions where it is more sensitive than RV measurements 
with 1 m s _1 precision. Meanwhile, in nonresonant regimes 
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Figure 2. The upper limit on the mass of an interior perturbing 
planet (upper panel) and exterior perturbing planet (lower panel) 
in the TrES-1 system from the 11-point analysis. The regions 
bounded by the contours correspond to a companion mass that is 
less than 100 (dotted), 10 (dashed), and 1 (solid) Earth masses. 

the RV method remains superior for a large portion of the 
parameter space. Additional transit timing data, particu- 
larly with precision that is comparable to the most precise 
of the given data, would lower the entire limit from TTV, 
provided that no secondary planet exists. Such data would 
render the TTV approach superior over a larger range of 
periods. 



4 DISCUSSION 

Ultimately a combined analysis of all available data, includ- 
ing studies of the stability of candidate systems, will provide 
the most robust and sensitive determination of the parame- 
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Figure 3. A comparison of the limits expressed in this paper for 
the 11-point analysis (solid curve) and the expected limits that are 
achieved with 11 randomly selected radial velocity measurements 
that have a precision of 5 m s — 1 (dotted curve) and 1 m s _1 
(dashed curve) for a secondary planet with an eccentricity of 0.05. 
The radial velocity curves are given by equation J5J using A\ 2 = 
9, N = 11, Q = 0, and Q = 2. Near the 4:3, 3:4, and the 
1:2 mean-motion resonances the TTV analysis of these data can 
place tighter constraints (about a factor of 3) on the mass of 
putative secondary planets than can the RV technique with 1 rn 
s precision. 

ters of each planetary system. None of the planetary systems 
that compose the limits in Figures 0] and 03 below 10M® 
are stable for more than 10 6 orbits; though stable orbits 
with comparable masses, periods, and eccentricities exist. 
An overall stability analysis to accompany these TTV analy- 
ses was prohibitively expensive. We estimate that for Figure 
|5|there were ~ 10 7 potential systems that were analysed and 
~ 10 systems were analysed for each curve in Figures Eland 
[^-requiring a total of ~ 10 s hours of processor time. General 
stability analyses are often better s uited to confirm o r stud y 
a particular ca ndidate system (e.g. lBarnes fc Ouinnl l)2004) : 
|ji et al.l (120051) ') than to constrain candidate systems in the 
manner that we have done in this section. We did conduct 
a stability analysis for the system shown in Figure Q and 
found it to be unstable within 10 6 orbits; no neighboring 
stable systems of comparable x 2 were found. 

It is unclear what fraction of probable planetary com- 
panions are excluded by our results. Aside from the fact 
that many of the known, mult i- planetary systems are in 
m ean-motion orbit a l resonance (e.g. GJ 876), recent works 
bvlThommesj (l2005h . |Papaloizou fc Szuszkiewiczl J2005I) . and 
IZhou et alJ 1120051) show that, under fairly general initial 
conditions, small planets are readily trapped in low-order, 
mean-motion resonances with a gas-giant as it interacts 
with a protoplanetary disk and migrates inward. These re- 
sults may imply that resonant systems are common among 
multiple-planet systems. If this is true, then results like 
Figure [5] actually exclude a much larger fraction of the 
probable orbits than one might infer from the projection 
of the excluded regions ont o the a/e plane. In addition, 
iMandell fc Sieurdssonl (120031) show that a large fraction of 
existing, terrestrial planets can survive the migration of a 
Jupiter-mass planet, though only a fraction of the survivors 
will be in resonance. 

The sensitivity of TTV to the mass of a perturbing 
planet renders it ideal for discovering and constraining the 
presence of additional planets in transiting systems like 
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Figure 4. This shows the perturber mass limit from a more re- 
fined TTV analysis near the 2:1 mean-motion resonances as a 
function of the period ratio of the perturber and TrES-lb. The 
upper panel is for an interior perturber and the lower panel is for 
an exterior perturber. The solid curve is for a perturber eccen- 
tricity of 0.00, the dashed curve is for an eccentricity of 0.02, and 
the dotted curve is for an eccentricity of 0.05. The long-dashed 
line shows a mass of 5Mg, the dashed line is for 2Mgj, and the 
dotted line is for lMj. For a perturber with eccentricity 0.02, the 
lowest point on the curve corresponds to a mass of 0.65Mg). 



Figure 5. This shows the perturber mass limit from a more re- 
fined TTV analysis near the 3:2 mean-motion resonances as a 
function of the period ratio of the perturber and TrES-lb. The 
legend for this plot is the same as in Figure HI 
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TrES-1. These studies can help determine the ubiquity of 
multiple planet systems and resonant systems-including the 
distribution of mass in those systems. Moreover, TTV anal- 
yses of several systems can play a role in identifying the im- 
portance of various planet-formation mechanisms. For ex- 
ample, the presence of close-in terrestrial planets favors a 
sequential-accretion model of planet formati on over a grav- 
itational instability model JZhou et al.ll2005tl . 

For the case of TrES-1, while a companion planet with 
a mass larger than the Earth is ruled out in a portion of 
the available parameter space, there remain large regions 
where additional planets could reside; this includes mutu- 
ally inclined orbits which were not studied in this work (we 
believe that these results are valid provided that the angle of 
mutual inclination i satisfies cosi ~ 1). The 11 timings anal- 
ysed in this work, excepting the questionable "12 th " point, 
deviate from a constant period by an amount that is difficult 
to interpret; having a \ 2 of about 2 per degree of freedom. 
Additional observations, with higher precision, will allow us 
to more fully constrain the system and to interpret the ex- 
isting transit timing variations. 
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